The model originally proposed by Rosenthal and Henderson (the RH model) and later elaborated by Liao and Frank connects several key factors in the CryoEM method using the following equation:^[26] 1 where N is number of particles (projections) that are needed to reach the resolution defined by the Fourier frequency k; σ² is the variance of the noise, is the Fourier intensity of the map in the resolution shell [k,k + Δk], B is the effective temperature factor, and C is a scaling constant.

The SNR (defined as the ratio between variances of signals and noises) at a resolution shell k is effectively represented by , where the numerator resembles the rotational average of the intensity (similar to the small/wide angle scattering intensity) and the denominator is the noise level. The second term on the right side of Eq. (1), k, describes the linear dependency of N(k) on the resolution shell k because the number of data points in each two-dimensional (2D) projection is at the order of k², while the desired number of data points in 3D increases with k³. The last term is the Gaussian falloff to account factors, such as structural fluctuations or misalignment during data processing. We carried out numerical simulations to examine the effects of each component. Based on the results, we revise the formula to provide an improved prediction formula on map resolutions.

The structure of GroEL (PDB ID: 1XCK,^[34]) was used as the model system in the numerical simulations (see Fig. 1). We present two theoretical formulations to investigate the factors that influence the map resolutions, the Gaussian noise (GN) model and the Thermal Fluctuation (TF) model, as shown in Fig. 1.

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. (color online) The two models used for CryoEM single particle imaging data simulations. (a) Gaussian noise model: noises were added to the simulation data to study the influence of noise on model resolution. (b) Thermal fluctuation model: a structure ensemble was first generated to mimic the structure heterogeneity. Each single particle projection was simulated based on a randomly picked structure from the ensemble. Three representative models for GroEL are shown in panel (b): the original structure is in gray; the yellow and blue structures were generated using normal mode perturbations.

The Gaussian noise (GN) model was proposed to describe the dependency of the model resolution on the number of projections and the noise variance. 2 where N(k) is the desired number of particles (projections) to collect confident signal levels at a frequency k; σ² is the variance of the Gaussian noise; and B is a scaling parameter.

The TF model describes how thermal fluctuations of the molecules influence the model resolution. RMSD (root-mean-square-deviation) is one quantity that measures the difference between the structures, and here we used the mean square of the pairwise RMSD of an structure ensemble to quantify the structure fluctuations and to mimic the Debye–Waller factor. The following formula is proposed to relate the thermal fluctuation levels, number of projection images and achievable map resolution: 3 where 〈 RMSD²〉 is the mean square of the RMSD values obtained by pairwise comparison within the structure ensemble; and k is the spatial frequency. In both Eqs. (2) and (3), A_GN, A_TF, B, and C are free parameters.

For the Gaussian noise model, SPIDER 22.03^[35] packages were used to generate the simulation data. The atomic model of GroEL was first converted to a density map (voxel size = (0.86 Å)³), and then projection images were simulated at orientations generated with the successive orthogonal rotation sample approach.^[36] The noises were incorporated according to the desired SNR following a Gaussian distribution after the contrast transfer function (CTF) for defocus of range from 1.0 μm to 3.0 μm was convoluted to the simulated projections. In this case, the CTF was not modulated by the envelope function. The image simulation process is summarized in Fig. 2(a). More specific procedure in mathematical phrase is as the following: (i) a 3D map is rendered by voxelization of an atomic model; (ii) generate projection images at orientations that uniformly distributed in the SO3 space; (iii) convert the projection images to Fourier space using FFT, then multiplies the CTF function for a random defocus level between 1.0 μm to 3.0 μm; (iv) apply inverse Fourier transform to get the CTF convoluted projection image; and (v) add desired Gaussian noise to get the final projection images.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) The workflow for single particle imaging data stimulation. (a) The protocol for the projection data simulation with defocused lenses and Gaussian noises from a single 3D structure. (b) The procedure for structure ensemble generation to model the heterogeneity of biomolecular conformations.

To simulate the heterogeneity in the structures, we first obtained a set of diverse structures to form a structure ensemble based on the GroEL structure. Without the loss of generality, the structural ensemble was generated using the normal mode perturbation approach. The ‘ProDy’ program based on an anisotropic elastic network model was used to compute the normal modes and the eigenvalue spectrum.^[37] Since the functional relevant motions are highly collective, three normal modes corresponding to the lowest frequencies were used to generate the perturbed structures. Specifically, the original structure (gray colored in Fig. 1(b)) was perturbed along the three normal modes (accounting to approximately 20.8% of total fluctuations), with deformation amplitudes ranging from 1 to 100, and an ensemble of 1000 structures was generated around the original structure, with a maximum RMSD of approximately 10 Å compared to the original structure. The RMSD values of generated structures with respect to the original structure were used to group the structures into 10 bins using 1 Å as the bin size. The structure ensembles were then compiled by drawing structures from the 10 bins using the following approach: the first group includes structures randomly drawn from the first bin (RMSD < 1 Å compared to the original structure). The second group is composed of structures randomly drawn from the first two bins, etc. This procedure ensures that the structure diversity levels are different in 10 groups. The structure preparation procedure is summarized in Fig. 2(b). The average structure deviation from the original structure and the mean square of the pairwise RMSD within each group, which is used to measure the structure diversity, are summarized in Table 1. This allows us to examine the dependence of the model resolution on the structure heterogeneity (or thermal fluctuation levels).

Table 1.

Table 1.

Table 1. The characteristics of the structure ensembles after incorporating the thermal fluctuations. In each group, there were 100 structures selected using the protocol described in the main text. In the table, 〈RMSD〉a: Average RMSD with respect to the original model. 〈RMSD2〉b: The mean square of the pairwise RMSD within the ensemble. . Group ID 〈RMSD〉a/Å 〈RMSD2〉b/Å2 1 0.79 0.90 2 1.24 1.65 3 1.62 3.16 4 2.08 5.86 5 2.52 10.78 6 3.02 17.35 7 3.54 23.03 8 4.02 27.88 9 4.53 32.51 10 4.99 35.80

Table 1. The characteristics of the structure ensembles after incorporating the thermal fluctuations. In each group, there were 100 structures selected using the protocol described in the main text. In the table, 〈RMSD〉a: Average RMSD with respect to the original model. 〈RMSD2〉b: The mean square of the pairwise RMSD within the ensemble. .

The projection image simulation procedure is essentially the same as previously described for the case of a single 3D structure, except that images forming a dataset were simulated based on randomly selected 3D structures from the corresponding group. Consequently, each simulated dataset has the following controlled parameters: the number of projections, SNR, and the structure heterogeneity due to thermal fluctuations.

The resolution determination was based on the Fourier shell correlation (FSC) implemented in SPIDER with the gold standard rule at the cutoff level of FSC = 0.143. The simulation data was split into two half subsets randomly, and each reconstruction was carried out using back projections with known orientations or using the iterative reconstruction methods based on Bayesian maximum likelihood algorithm implemented in the Relion 1.4 package.^[4] In the cases that the standard reconstruction procedures were carried out to build the electron density maps, each dataset was processed five times and the best resolutions were used in the final analysis.

The resolutions of the reconstructed maps were determined for various combinations of factors that were considered in this work. The value of each factor was systematically scanned in practical ranges so that the quantitative relationships could be studied using a parameter fitting to theoretical formula. The free parameters used in the GN and TF models were obtained by the nonlinear Least Squares (Curve Fitting) module in MATLAB.

A statistical survey was carried out on the structures determined using CryoEM deposited in the EMDB database (http://www.ebi.ac.uk/pdbe/). The retrieved information includes the molecular weight, detector type, resolution determination criteria, molecular symmetry, number of projections and the microscope operating voltages.

We focused on the models that fulfill the following requirements: (I) Models that were deposited between 2017/01/01 and 2018/05/17; (II) Micrographs were recorded using direct electron detection technology; (III) The reported resolution is determined with gold standard rule at FSC = 0.143; (IV) Models without higher symmetry; (V) Molecule weight of all models are ranged from 0.5 MDa to 1.0 MDa. As a result, 86 models from EMDB were selected to for the map resolution statistics.

The RH model describes the dependency of map resolution on the number of experimental images, and the logarithm trend of the reconstructed model resolution as a function of the number of measured projections is attributed to the b-factors (due to sample particle drifting, misalignment, numerical interpolation, etc.^[26] Surprisingly, the logarithm trend was observed for the simulation data without explicitly applying the envelope function e^Bk2/2 during the projection simulations (see Fig. 3(a)). Furthermore, the utilization of the known orientation as available information eliminated alignment errors that occurred during orientation recovery process. Therefore, the errors described with the b-factor envelope function in RH model are not the only cause of the logarithm trend. The numerical interpolation during structure reconstruction could not be avoided, and this might be one source of the logarithm trend. Another possibility is that the logarithm trend is due to uneven information embedding in single particle imaging datasets, where the low-resolution information is redundant, while the high-resolution information barely reaches the signal-to-noise threshold.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. (color online) Relation between model resolution and data quantity/quality. (a) The model resolution depends on the number of measurements. Two sets of data are plotted with signal to noise ratios (SNR) of 0.2 and 0.4. The lines were obtained by fitting the data points using the Gaussian noise model. R-square values were 0.9946 and 0.9865, respectively. (b) The noise term was combined with the number of projections. The two sets of data in Fig. 3(a) were fitted with a single set of parameters by defining the effective number of projections. R-square of the line is 0.9782. (c) and (d)The same plots as in panels (a) and (b), except that the maps were reconstructed using standard procedures in Relion without using known orientation information. Due to the heavy computing demands for standard reconstructions, a subset of data were investigated in the case of panels (c) and (d), and the sampled region in panel (d) corresponds to the blue shaded area in panel (b).

In the 2D cases, n measurements of the same image will boost the SNR (defined as ratio of signal variance and noise variance) n times if the noise types and levels are the same for all measurements.^[30] In the 3D map reconstruction from 2D projection images, we observed the same relationship (see Fig. 3). To simplify the GN model, a new parameter, the ‘effective number of projections’, N_e, was defined as the product of the projection number and the average SNR of the individual projection image (see Eq. (4)). The noise term is absorbed into this new parameter N_e, and equation (2) becomes: 4 This relation can be verified with the simulation results summarized in Fig. 3.

In Fig. 3(a), the SNR values were treated as a separate parameter, independently from N (the number of measurements) and the two sets of data with different SNR values were fitted to two equations (the two black curves). In Fig. 3(b), the data points were merged to a single line using the effective number of projections (N_e). Subsequently, a single equation is adequate to describe the relationship between the resolution and the number of ‘effective’ projections. Using the nonlinear curve fitting algorithm, the values of coefficient A_GN and B are 1298 and 46 (Ų), respectively, for the GroEL simulation data.

The map reconstructions using back projection method are only possible if the image orientations were known, so it describes an ideal situation. In practice, the orientations needs to be assigned using iterative methods, which converges to the model that best satisfies the constraints of the whole dataset. Here, the auto-refine function implemented in Relion was applied for the model reconstruction. The gold standard rule was also used for resolution determination. To reduce the influence from randomness, the best resolution from five independent runs were reported as the map resolution for final statistics. The results are summarized in Figs. 3(c) and 3(d). The orientation assignment error indeed affects the quality of the reconstructed maps, indicated by the lower resolutions compared to that in the cases from back projection reconstructions. Nevertheless, after combining the number of projection images and the SNR, the overall trend is still reasonably described using the RH model, suggesting that the effective number of particles is valid in the simulated cases. However, the relation becomes invalid if the SNR is too low for accurate orientation assignment in extreme cases, which were not considered in this study.

Biological macromolecules exist in a thermal environment, and the structure fluctuates around the native states. In many cases, due to the functionality, molecules exist in several meta-stable conformations.^[38–40] Using the normal mode perturbation approach, we attempted to simulate intrinsic motions and study their influence on the achievable structure resolutions at various conformation heterogeneity levels. The TF model described by Eq. (3) captures the relationship between map resolution and the number of projections in the presence of thermal fluctuations (i.e., structure heterogeneity), using an analytical formula similar to the case of the temperature factor in x-ray crystallography.^[41] Compared to the Gaussian noise model, the TF model predicts a different behavior in the resolution dependency on the number of measurements: the resolution gets worse faster for structure ensemble with larger thermal fluctuations.

In practice, both experimental noise and molecular thermal fluctuation have impacts on the CryoEM single particle experimental data. Therefore, it is necessary to derive a model that combines the GN and TF models. Intuitively, the following formulation is devised by treating the Gaussian noise and thermal fluctuation as independent factors that affect the structure resolutions: 5 Note that the effective number of particles is used in Eq. (5). To get a sensible fitting to the simulation data, we divided the TF model to two regimes using 〈RMSD²〉 = 10 Ų as a threshold. As a result, two sets of parameters were obtained: for 〈RMSD² 〉 > 10 Ų, A_TF = 1702.0, B_TF = 58.33 Ų, and C_TF = 0.75; for 〈RMSD² 〉 ⩽ 10 Ų, A_TF = 2465, B_TF = 41.06 Ų, and C_TF = 2.88. Without the thermal fluctuation term C 〈RMSD²〉 in Eq. (5), we will need five sets of parameters to fit five datasets in Fig. 4.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. (color online) The influence of thermal fluctuation to the structure resolution. (a) The relationship between the resolution and number of projections at various thermal fluctuation levels. The thermal fluctuation is quantified using the average square root-mean-square-derivation 〈RMSD2 〉 of the structure ensemble. All the projections were simulated with Gaussian noise with an SNR of 1.00. The R-square values of the fittings to the GN model are approximately 0.96. (b) The five datasets were fitted using a single set of parameters with an R-square of 0.98 using Eq. (5).

To compare our theory and simulation results with experimental data and cross-validate the conclusions, we conducted a systematic survey on the structures determined using CryoEM single particle imaging technology, which are deposited in the EMDB database. It is worth noting that only the subset of structures that met the criteria described in the Methods section were used for the statistics. The resolution distribution nicely resembles the relationship described using the noise model and the thermal fluctuation model. Using these parameters (A_GN, A_TF, B, C in Eqs. (4) and (5)), we can draw theoretical guidelines to estimate the required number of projections and obtain the desired resolutions; see Fig. 5. Because the original model 1XCK weighed approximately 0.81 MDa, the data for molecule weight ranged from 0.5 MDa to 1 MDa were compared with the analytical curves of the two models.

Fig. 5.

Figure Option ViewDownloadNew Window

Fig. 5. (color online) The distribution of model resolutions as a function of particle numbers for experimentally determined structures. The theoretical lines indicate the resolution limits for molecules with molecular weights between 0.5 MDa and 1.0 MDa. The curves were drawn with SNR = 0.01 and 0.18 for upper and lower bounds predicted using GN model. The upper bound were predicted using TF model with 〈RMSD2〉 = 50 Å2 and SNR = 0.01, the corresponding lower bound were predicted with 〈RMSD2 〉 = 0 Å2 and SNR = 0.18.

As shown in Fig. 5, the experimental model resolutions were bounded by the theoretical curves. In the bounds predicted by GN and TF models, the parameters obtained from GroEL data fitting were used. For GN model, the SNR levels were estimated to be 0.01 and 0.18 for upper and lower bounds. For the TF model, the same SNR levels were used, except that the upper bound has a thermal fluctuation term with 〈RMSD² 〉 = 50 Ų. The lower bounds for both GN and TF models are very similar; however, it is clear that the TF model provide a better estimate of the upper bound of resolution for experimental data. Another interesting observation is that the resolution hardly broke the barrier of 3 Å, although the number of projections covers a broad range from 2 × 10⁴ to 2 × 10⁵. This is likely due to the thermal fluctuations or the coexistence of multiple conformations of the molecule.